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Abstract: 

We  are  studying  here  the  numerical  simulation  of  high  Reynolds  number  internal,  2-D  flow  in  a 
convergent-divergent  nozzle,  for  a  compressible,  viscous  fluid.An  implicit  finite-difference  scheme  is 
used  to  solve  the  parabolic  approximation  of  the  Navier-Stokes  equations,  so  that  the  time  steps  are 
not  severe'y  limited  by  the  small  grid  sizes  needed  for  the  computation  of  the  vicous  effects.  After 
the  resolution  of  this  problem  on  a  serial  computer,  we  describe  the  first  steps  of  the  parallelization 
of  this  problem  on  a  Hypercube  (parallel  version  of  the  ADI  method  (7,10)). 
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1.  Introduction: 

The  resolution  of  the  fluid  flow  in  a  nozzle  is  a  problem  which  has  been  studied  by  several 
authors  on  sequential  [3,4,9,15]  and  vector  [14]  computers.  At  the  present  time,  no  complete  study 
of  a  fluid  dynamic  problem  has  been  made  on  a  distributed  memory  parallel  machine  such  as  the 
iPSC  Intel  Hypercube,  and  we  propose  here  to  study  the  parallelization  of  this  physical  problem 
and  its  corresponding  speedup  on  such  a  machine. 

A  general  technique  of  resolution  of  the  Navier-Stokes  equations  [12,13]  mainly  based  on  the 
combination  of  a  general  curvilinear  coordinate  transformation  and  an  implicit  method  has  given 
some  good  results.  It  has  already  been  experimented  on  the  ILLIAC  IV  [8]  and  will  be  used  here. 

The  physical  domain  (the  domain  inside  the  nozzle)  is  first  transformed  into  a  rectangular 
domain  by  the  resolution  of  a  pseudo-elliptic  equation  with  Dirichlet  conditions  at  the  boundaries 
(Thompson  et  al  method  [16]). 

The  Navier-Stokes  equations  are  then  solved  in  this  domain  with  the  Beam- Warming  implicit 
method  [1],  which  leads  to  the  computation  of  two  sets  of  linear  systems  (ADI  method). 

The  modelization  of  the  physical  problem  is  descibed  in  section  II,  and  the  numerical  generation 
grid  in  section  III.  The  boundary  and  initial  conditions  are  traited  in  section  IV,  and  we  deal 
with  the  turbulence  modelization  in  section  V.  The  numerical  algorithm  is  explained  in  section 
VI,  the  results  and  their  discussion  appear  in  section  VII.  Finally  in  section  VIII,  we  discuss  the 
parallelization  of  this  problem  on  a  Hypercube. 

2.  Modelization: 

2.1.  Navier-Stokes  equations: 

The  motion  inside  the  nozzle  (see  figure  1)  is  governed  by  the  unsteady  averaged  Navier-Stokes 
equations  in  the  inertial  cartesian  coordinates  (x,y,t).  These  equations  written  in  the  non-dimension 
form  and  with  the  strong  conservation  law  form  [12,17]  are  given  by: 
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and: 

p:  density. 

u,  v:  cartesian  velocity  components. 
en:  total  energy. 
p:  pressure. 

a:  sound  speed  (o  =  y/'i RT). 

T :  temperature. 

7:  ratio  of  specific  heats  (dry  air:  7  =  1.4). 

/*:  dynamic  viscosity. 

A:  taken  with  the  Stokes  hypothesis  A  =  -(2/3)/*. 
k:  coefficient  of  thermal  conductivity. 

R:  gas  constant. 

U:  modulus  of  the  velocity  (17  =  %/tt2  4-  uJ). 

D:  specific  length. 
ep:  specific  heat. 

Re:  Reynolds  number  (Re  =  (UDp)/p). 

Pr:  Prandtl  number  (Pr  =  ( pevjk ). 

M:  Mach  number  (M  =  Ufa). 

The  equations  (1)  and  (2)  are  completed  by  the  equation  defining  the  pressure: 
p  =  (7  “  !)[«*  ~  0.5p(uJ  4-  w2)] 

As  our  nozzle  is  symmetric  about  the  central  axis  y=0,  we  can  restrict  our  study  to  the 
side  of  the  nozzle. 


Figure  1:  45°  -  15°  Convergent- Divergent  Nozzle 


2.2.  Transformed  Navier-Stokes  equations: 

2.2.1.  General  form: 

The  resolution  of  equations  (1)  in  the  region  (R),  closed  by  the  boundary  (S)  (see  figure  (1)) 
can  be  greatly  simplified,  by  transforming  these  equations  (l)  to  a  new  body-fitted  curvilinear 
coordinate  system  [15,16]. 

More  over,  if  we  take  for  the  dependent  variables,  the  Cartesian  velocity  components,  the  basis 
equations  (1)  can  be  transformed  under  the  independent-time  invertible  transformation: 

(£  =  £(*>y)  (*  =  *(£,»?) 

<»?  =  ^(i,y)  4=>  y  =  y((,i?)  (5) 
l  r  =  t  W  =  r 

to  an  arbitrary  curvilinear  coordinate  space  (£,  »?,r)  while  maintained  the  strong  conservation 
law  form  of  eq  (l)  [8,12,17].  The  next  chapter  will  present  the  numerical  generation  of  the  curvi¬ 
linear  space  (£, refered  as  the  computational  domain. 

The  resulting  equations  are  given  by: 

drq  +  d(£  +  dnf  =  Re~1[d((J~1(Zx9i  +  ^2))  +  +  ^<72))]  (6) 

In  our  particular  physical  case,  due  to  the  geometric  properties  of  the  nozzle,  we  can  restrain 
the  transformation  (£,r),r)  to  [14]: 


£  =  £(*) 
t|  =  n(i,y)  (?) 

r  =  t 


The  equations  are  then  rewritten  as  follows: 


dTq  +  dte  +  dr,f  =  Re  1[di{J~1{^xgi))  +  d^J  1(r)xg1  +  rjyg2))\ 


where,  the  vectors  q,  e  and  /  are  expressed  by  the  relations: 


q  =  J~lq 


(8) 


(9) 


By  the  definition  of  the  velocities: 


e  =  J~l{Zze) 

f  =  J~1(Vze+  r]vf) 
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called  the  contravariant  velocity  components,  and  corresponding  to  the  decomposition  of  the 
vector  velocity  along  the  £  and  r\  curvilinear  coordinates,  we  get  the  following  expressions: 
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J  is  the  transformation  Jacobian: 

J  =  tzVv  =  l/(*€Vr >)  (11) 

The  cartesian  derivatives  such  as  uz  are  expanded  in  the  (£,  rj)  space  via  chain-rule  relation: 

t*z  =  +  W*  (12) 


And,  the  metrics  formed  from  chain-rule  extension  of  if, y^,  y,,  are  given  by  the  fol¬ 

lowing  relations: 


Zz  =  Jyr, 

Tjz  =  -Jyt  (13) 

f?„  =  Jzf 


2.2.2.  Simplified  form:  the  parabolic  equations. 

Classicaly  for  high  Reynolds  number  flows,  we  solve  the  viscous  terms  only  near  the  rigid 
boundaries.  We  also  make  the  hypothesis  of  the  parabolic  approximation  that  the  viscous  terms  in 
(  (along  the  body)  are  neglected  and  only  the  viscous  terms  in  q  are  retained. 

The  equations  (8)  are  then  equal  to  [8,12,15]: 


dTq  -I-  +  dnf  =  Re  ldng 


(14) 


where: 
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Note  that  unlike  boundary-layer  theory,  the  pressure  p  can  vary  through  the  viscous  layer,  and 
all  the  inertial  terms  of  the  normal  momentum  equation  are  retained. 
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3.  Numerical  grid  generation: 

By  the  means  of  the  generation  of  (£,  >?,r),  the  physical  domain  is  transformed  into  a  rect¬ 
angular  domain  with  a  square  grid  (A£  =  A rj  =  1).  All  the  computations  are  then  done  in  this 
rectangular  domain  (the  computational  domain)  [16],  using  more  simple  finite-difference  operators 
for  the  different  derivatives  of  the  Navier-Stokes  equations. 

Note  that  the  other  advantage  of  this  transformation  is  that  the  boundary  surfaces  (here  the  wall  of 
the  nozzle)  are  mapped  into  rectangular  sufaces,  so  that  the  boundary  conditions  can  be  computed 
more  easily  and  more  accurately. 

Also,  this  transformation  allows  grid  point  clustering  near  the  walls,  such  that  the  viscous  effects 
are  well  computed  (see  figures  2  —  3). 

The  numerical  method  of  generation  of  the  grid  is  given  here  by  the  Thompson  scheme.  In  this 
method,  the  grid  in  the  physical  domain  is  determined  by  the  solution  of  a  Laplace  or  a  Poisson 
equation  [12,16]. 

d2£/dx2  +  d2£/dy2  =  0  or  /(£,»;)  (I6.a) 

d2r)/dx2  +  d2rj/dy2  =  0  or  /($,»?)  (16.6) 

where  the  values  of  (  and  rj  are  arbitrarily  fixed  on  the  boundaries. 

Pratically,  the  equations  (16)  are  transformed  into  the  computational  domain  and  the  generation 
of  the  grid  is  then  given  by  the  resolution  of  the  pseudo  -elliptic  equations: 

ax£€  -  20z€„  +  - ixtl =  0  or  -J2{txxy(  4-  /(£,  fj)y„)  (17. a) 

ay(e  -  2pyin  +  =  0  or  -J2{£xx y€  +  /($, rj)y„)  (17.6) 

In  our  case,  we  have  arbitrary  chosen  i  =  z(£)  as  an  exponential  law  centered  in  the  throat. 
The  coefficients  of  the  equation  (17.6)  to  solve  are  [14]: 

«  =  !/*>  &  =  J/CVn.  7  =  +  y|,  and  £xx  =  -x££/z|  (18) 

The  boundary  values  of  y  are  known  and  correspond  to  the  values  of  the  desired  mesh  point 
on  the  boundaries  of  the  physical  domain. 

The  function  is  determined  at  every  mesh  point  (£ ,rj ),  and  its  function  is  to  control  the 

spacing  between  the  grid  points  in  the  interior  of  the  physical  domain.  Here,  we  have  chosen  the 
function  given  by  Strikwerda  in  [14]. 

The  following  figures  (2  —  3)  show  the  grid  lines  £  =  £Const,  and  rj  =  t]con»t  in  the  physical  domain 
for  the  Laplace  and  Poisson  equations. 


Grid  lines  (£,  rj)  in  the  physical  domain  (  for 
ation  ). 


4.  Boundary  conditions-  Initial  condition: 

The  equations  (14)  have  to  be  solved  in  the  computational  domain  defined  before,  subject  to 
different  conditions  on  the  boundaries. 

The  curvilinear  coordinate  system  (£,rj)  is  defined  such  that  the  inflow  and  outflow  boundaries  are 
two  (-coordinate  lines,  and  the  wall  and  the  symmetric  axis  are  two  ^-coordinate  lines. 

4.1.  Wall  boundary  condition: 

At  the  nozzle  wall,  we  have  for  the  velocity  the  no-slip  condition: 

u  =  v  =  0  (19. a) 

and,  for  the  temperature  we  can  take  the  adiabatic  condition: 

n  >  VT  =  0  (19.6) 

with:  n  the  vector  normal  to  the  surface.  In  our  case,  for  a  surface  q  =  r)contt,  we  have: 

n=  Vq  (19. c) 

and  equation  (19.6)  becomes: 

V»7.Vr  =  0  (19. d) 

Under  the  parabolic  approximation  is  neglected  and  we  have  [15]: 

(Vq  •  Vq)r„  =  0  (19. e) 

4.2.  Symmetric  axis  y=0: 

Because  of  the  symmetry  of  the  nozzle,  we  consider  the  axis  y=0  as  one  of  the  boundaries  of 
this  problem  (see  figure  (1)). 

In  order  to  limit  the  overloading  of  the  code,  we  just  limit  ourselves  to  a  condition  on  the  vertical 
component  of  the  velocity  v  : 

v  =  0  (20) 

as  it  is  an  odd  function  of  y. 

4.3.  Upstream  inflow  boundary: 

The  pressure  p  and  the  temperature  T  are  constant  during  all  the  computation,  and  are  equal 
to  their  stagnation  correspond  ants  i-e: 


The  components  u  and  v  of  the  velocity  are  given  by  the  relation: 

u  =  t/tan0(y)  (21-c) 

where  6(y)  is  a  function  given  by  Holder  et  al  [5]. 

4.4.  Downstream  outflow  boundary: 

For  a  supersonic  flow  (M  >  1),  no  boundary  conditions  are  required.  The  variables  are  deter¬ 
mined  by  extrapolation  from  the  interior  points. 

4.5.  Initial  condition: 

All  the  variables  at  the  beginning  of  the  computation  were  computed  under  the  1-D  approxi¬ 
mation  method  (isentropic  flow,  perfect  fluid). 

In  such  an  approach,  the  values  of  the  variables  at  a  specified  section  of  the  nozzle  are  given  by 
the  geometric  ratio  of  the  aera  of  the  section  over  the  throat  section,  and  by  the  chosen  stagnation 
conditions  [6,11]. 

In  particular,  we  have  taken  here: 

p0  =  6.2  10 6  N/m2 
To  =  300 nr 

And  the  Reynolds  number  for  this  initial  condition  was  then  equal  to:  Re  =  1.5  106 


5.  Turbulence  Modelling: 

For  Re  >>  1,  the  domain  is  divided  into  two  main  regions:  the  boundary- layer  and  the  core. 
The  modelization  of  the  turbulence  has  to  separate  these  parts: 

5.1.  Internal  and  external  wall  boundaries  layers: 

The  boundary  layer  flow  along  the  wall  is  modelized  by  using  a  two-layers  eddy  viscosity  model. 
Under  this  approach,  the  turbulent  stresses  u  in  the  boundary-layer  are  modeled  in  terms  of  the 
eddy  viscosity  fit  by: 

rt  =  fitdU/dy  (22) 

where  U  is  the  velocity  component  parallel  to  the  wall,  and  y  is  the  coordinate  normal  to  and 
measured  from  the  wall  [2,15]. 

For  the  inner  layer,  the  eddy  viscosity  fit  is  given  by  the  Van  Driest  formulation  with  Cebeci 
damping: 


fit  =  0.016py*[l  -  exp(y+  /A+))2dU/dy 


(23) 


where  the  damping  constant  A+  is  function  of  the  pressure  gradient  parameter  P+  as: 

A+  =  26(1  +  11.8P+)-1/2  (24) 

where: 

P+  =  n/(p2u*).dp/da  (25) 

with  the  “friction”  velocity  given  by: 

tif  =  y/  Tw/ Ptti  (26) 


where  the  index  w  means  that  the  values  are  taken  at  the  wall. 

For  the  outer  layer,  the  eddy  viscosity  is  computed  from  the  Clauser  approximation: 


Hr  =  0.0168p 


[\ue-U)dy 

Jo 


(27) 


where: 

Ue  :  core  flow  velocity  at  the  edge  of  the  boundary  layer. 

6  :  thickness  of  the  boundary  layer. 

Pratically,  the  boundary-layer  thickness  8  is  defined  here  as  the  distance  from  the  wall  so  that 
the  velocity  U  approches  the  corresponding  value  of  the  freestream. 

5.2.  Core  flow: 

The  effect  of  the  turbulence  is  also  approximated  in  this  part  from  the  modelization  of  an  eddy 
viscosity,  constant  for  all  the  core  flow  and  deducted  from  the  jet  theory  [15]: 
we  have  here: 


fit  =  0.0256bpUo 


(28) 


with: 

U0\  the  average  inlet  velocity. 
b:  the  nozzle  inlet  radius. 

Finally,  a  constant  turbulent  Prandtl  number  of  Pr*  =  0.9  is  used  in  the  energy  equation, 
which  leads  to  the  following  definition  of  the  turbulent  eddy  conductivity  kt: 


Prt  =  fitcp/kt  =  0.9 


(29) 


6.  Numerical  Method: 

We  solve  the  equations  (14)  in  the  computational  domain  (£,»?)  subject  to  the  boundary 
renditions  seen  before,  with  the  implicit  delta-form,  approximate-factorization ,  Beam- Warming 
algorithm. 

An  implicit  numerical  method  is  employed  here  in  order  to  avoid  the  severely  restrictive  stabil¬ 
ity  conditions  of  an  explicit  method,  when  small  grid  spacings  are  used.  Such  a  situation  is  needed 
near  the  wall  for  an  accurate  computation  of  the  viscous  effects. 

In  the  basic  Beam- Warming  algorithm,  the  spatial  derivative  terms  in  the  Navier-Stokes  equa¬ 
tions  are  approximated  by  standard  second-order  accurate  central-difference  operators,  and  the 
implicit  5-method  of  Richtmyer  and  Morton  is  taken  for  the  time  differencing  [8] . 

The  computation  of  the  boundary  points  can  be  computed  directly  by  the  numerical  resolution 
of  the  equations  (see  [15])  or  by  extrapolation  from  the  interior  points  at  the  end  of  each  time  step 
(see  [12]).  The  second  case  degrades  the  time  accuracy  on  the  boundary  to  a  zero-order  but  gives 
a  more  simple  scheme. 

The  method  proposed  here  is  to  mix  these  two  approaches:  the  variables  at  the  inflow  boundary 
are  updated  by  extrapolation  (due  to  its  particular  conditions),  and  for  the  other  boundaries,  we 
use  a  direct  numerical  resolution.  This  approach  will  give  us  the  correct  steady  state  solution  but 
will  need  more  time  steps. 


6.1.  Time  Differencing: 

For  both  the  interior  and  boundary  points,  we  define  the  same  time  differencing.  The  equations 
of  Navier-Stokes  in  their  final  form  (14)  are  rewritten  here  as: 

dTq  =  ~[d(e  +  d„f  -  Re~ldng }  =  f  (30) 

Using  n  for  the  time  index  and  h  for  the  time  step,  we  apply  the  Richtmyer  and  Morton 
method,  and  we  obtain  [8]: 


qn+l  =  qn  +  *((1  "  9)rn+1  +  0fn ]  (31) 

This  method  is  first-order  accurate  in  time  for  9  =  0  (Implicit  Euler  method),  and  is  second- 
order  accurate  in  time  for  9  =  1/2. 

Since  we  seek  only  the  asymptotic  steady  state  solution,  we  can  employ  the  first-order  accurate  in 
time  method.  The  accuracy  of  the  solution  is  given  by  the  spatial  difference  operators  [1]. 

So  we  have: 


qn+ 1  =  qn  +  hfn+1 


(32) 


In  order  to  define  the  non-linear  term  ?„+ j,  we  must  locally  linearized  the  terms  e,  /  and  g  in 
terms  of  q.  This  is  done  by  using  the  Taylor  series  expansions: 


en+1  =  e"  +  £"(?n+1  -  qn )  +  0{h2). 


fn+ 1  =  fn  +  Fn(qn+1  -  qn )  +  0(h2). 


(33) 


ES 


$n+1  =  gn  +  Gn(qn+1  -  qn)  +  0{h2). 


where  E,  F  and  G  are  the  flux  Jacobian  matrices: 


E  =  de/dq,  F  =  df/dq  ,  G  =  dg/dq 


(34) 


defined  more  precisely  as  following: 


E  =  [£y  =  dti/dqj,  F  =  [Fitj]  =  dh/dq,  ,  G  =  [Gu]  =  dg{/dq}  (35) 

By  definition,  the  flux  vectors  e  and  /  are  both  linear  combinations  of  e  and  /  (see  eq(9))  i.e: 

f  =  J~l{r)xe  +  r}vf)  (36) 

so  the  inviscid  flux  Jacobian  matrices  E  and  F  are  expressed  as  follows: 


E  =  de/dq  =  £xde/dq 
F  =  df/dq  =  r\xde/dq  +  r]vdf/dq 


(37) 


The  E  or  F  matrices  are  then  given  by  (8,12]: 

E,F  = 

0  Kx  K2 

Kit2  -  u*  e-Kifr-  2)u  K2u  -  (7  -  1  )Kiv 

K2t2-v8  Kxv-K2{t-  1)u  8-K2( 7  -  2)v 

.6{2<t>2  -  q(en/p))  [Kx[~i(en/p)  -  <t>2\  -  (7  ~  1)«0|  [#2(7 («»»/p)  -  4>2)  -  (7  -  1  )v6] 

(38) 


^i(7-l) 
#2(7  -  1) 

7* 


where: 


=  0.5(7  -  l)(u2  +  v2) 


=  K\\i  +  K2v 


with  the  following  definition  of  the  coefficients  Kx  and  K2\ 


for  E:  K 1  =  £*,  #2=0 


for  F:  K\  =  t]x,  K2  =  rjy 


(39) 

(40.a) 

(40.6) 
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And  for  the  viscous  flux  Jacobian  term,  we  have: 


where: 


and: 


(  0 

0 

0 

0 

nn  —  t~1 

021 

<*iM1/p) 

a^5„(l/p) 

0 

\j  — 

931 

0125,,  (1/p) 

ctsdr,(l /p) 

0 

\9* 1 

9*7 

9*3 

“45„(l/p) 

=  - 

/  dn(u/p) 
Q  Wv/p) 

^ ,  with  Q 

_  /  011  012 

\  012  013 

oti  =  /*[(4/3)>j*  +  »?’] 
a2  as  (n/Z)r)xT)v 


(41) 


(42) 


(43) 


<*3  =  M[»7 1  +  (4/3)i7j] 
Q!4  =  ~tkPr'l{r)l  +  t)l) 


By  means  of  the  local  linearizations  (see  eqs  (33))  the  equation  (32)  is  then  rewritten: 


[7  +  h(d^En  +  d„Fn  -  Re~1drtGn)] Af*  =  hfn 


(44) 


with: 


Ag"  =  gn+I  -  qn  (45) 

where  the  linear  operator  notation  is  to  be  interpreted  as  following: 

(d{£7*)A?n  =  di(En&qn)  (46) 

The  left  hand  side  of  equation  (44)  can  be  factorized  into  a  product  of  two  one-dimensional 
operators,  with  the  same  order  of  accuracy  [1]. 

This  results  in  the  factorized  form  of  equation  (44): 

([/  +  +  h(d„Fn  -  tfe-^G")]) Aqn  =  hfn.  (47) 

which  is  said  to  be  the  “delta-form”,  because  the  left  hand  side  contains  the  factor  A qn. 

The  r.h.s  of  equation  (44)  is  defined  as: 

r"  =  -fd€e"  +  d„fn  -  Re-ldr,gn\.  (48) 

Finally,  the  vector  solution  qn+1  i3  given  by  the  following  ADI  sequence: 


[I  +  hd€En]Aqtn  =  hfn 
[I  +  h(drtFn  -  Re-1dnGn)]Aqn  =  A q*n 


(49.a) 


(49.6) 

qn+l  =qn  +  A  qn  (49  .c) 

For  convenience,  we  have  omitted  the  spatial  subscript  notations,  all  the  terms  were  taken  here 
at  the  point  (*,j)  of  the  computational  domain. 

But  from  now  on,  we  will  use  the  spatial  subscripts  again. 

Note  that  by  the  choice  of  a  central-difference  scheme  for  the  space  derivatives,  each  step  of  the 
ADI  sequence  (49)  involves  the  solution  of  a  linear  system  of  equations  having  a  block-tridiagonal 
coefficient  matrix. 

6.2.  Space  Differencing: 

6.2.1.  Interior  points: 

The  resolution  of  the  ADI  sequence  (49. a,  49.6, 49. c)  is  completed  by  the  choice  of  the  finite 
difference  operators  S  for  the  spatial  derivatives  d $  and  dn. 

We  use  here  to  approximate  the  convective  derivatives  the  standard  central-difference  second 
order  accurate  operator.  For  example  we  have  for  the  first  derivative  of  equation  (49. a): 


SiEij  =  {dEtj/dt)  =  ( Ei+1J  -  Ei-ij)/ 2A£  (50) 

where  by  definition  of  the  grid  in  the  computational  domain:  A£  =  Ar)  =  1.  From  now  on, 
the  increments  A£  and  Arj  will  be  replaced  by  1  in  all  the  spatial  difference  formulas. 

First  we  get  the  expressions: 


fyEij  —  (Ei+ij  ~  Ei-ij)/2 

SjFij  =  (Fij+i  -  Fij-x)/ 2.  (51) 

From  the  general  form  of  the  viscous  term  written  as  follows: 

d[(oti,jd0ij/dr))\/dri  (52) 

we  see  that  we  have  to  approximate  the  viscous  derivative  in  a  more  complicated  way. 

If,  we  define  the  following  central-difference  operator  for  a  general  viscous  term  h: 


(dh,-j/5r?)  «  hij+i/2  hij— 1/2 


(53) 


its  application  on  the  viscous  term  defined  in  equation  (52)  gives  the  generalized  three-points 
central-difference  second-order  accurate  operator: 


d[(<*ij90ij/dr. l)]/dri  = 

l/2[(ai,;+i  +  —  (Qtj+i  +  2a,  j  -(-  a,  y 


-l  )Pi,j  +  (<*«',;  +  a.,j-i)A-,y-i] 
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(54) 


corresponding  to  the  following  notation  :  S'^gij  and  SjGij 

From  these  definitions,  the  ADI  sequence  is  rewritten  for  the  interior  points  as  follows: 


[/+A*js3yA«5  =  Mft 

(55.a) 

[/  +  A(«f -  R.-1«;g'V)]A£V  = 

(55.6) 

-V 
**•  + 

II 

•Q> 

C3 

+ 

> 

a 

(55. c) 

(56) 

6.2.2.  Boundary  grid  points: 

At  the  inflow  boundary,  both  the  pressure  p,  the  temperature  T,  and  the  density  p  remain 
constant  during  all  the  computation.  To  take  advantage  of  this  characteristic,  we  assume  that  for 
all  the  points  of  this  boundary,  the  increment  A qn  is  equal  to  0,  during  the  computation  of  the 
ADI  sequence.  Then  u,  v,  and  en  are  updated  by  extrapolation  from  the  interior  points  at  the  end 
of  each  time  step. 

For  the  other  boundaries,  we  used  a  method  of  resolution  quasi-similar  to  that  of  the  interior 
points.  In  this  case,  we  have  to  adapt  certain  derivative  operators  in  order  to  use  only  the  points 
of  the  computational  domain. 

We  will  employ  here  the  forward  difference  operator  A,  and  the  backward  difference  operator 
V  where: 


ij  —  *»+lj  — 

~  *»,;+!  —  *».i  (57) 

T7  AAA 

V«* ij  ~  tiJ  ~ 

This  approach  has  been  described  .  i  detail  in  [15],  and  is  presented  below  in  a  synthetic  way: 

6.2.2.I.  Outflow  boundary  (for  supersonic  condition): 

Only  the  first  step  of  the  ADI  sequence  (49. a)  has  to  be  modified.  The  spatial  derivatives 
along  the  (-axis  are  approximated  by  the  backward  difference: 

=  Eij  -  Ei-ij  (58) 

Due  to  the  link  between  the  finite-difference  and  finite  volume  methods,  we  associate  each  point 
of  the  outflow  boundary  with  a  half-cell  situated  A(/2  before  this  boundary.  The  time  derivative 
is  applied  at  the  centroid  of  this  half-cell  (placed  A(/4  from  the  outflow  boundary).  The  first  step 
of  the  ADI  sequence  is  then  rewritten  [15]: 


[/  +  fcV,(^  -  aI)]Aq%  =  hr *  (59) 

where  in  order  to  keep  the  second  order  accuracy  in  space,  the  coefficient  a  is  chosen  as  follows: 

a  =  1/(4  h)  (60) 


with  the  r.h.s  of  eq  (59)  given  by: 


fa  =  +  %  fa  ~ 


(61) 


6.2.3.  Wall  boundary: 

The  method  here  is  to  approximate  the  spatial  derivative  along  the  rj-axis  by  the  backward 
difference  operator  V  for  the  convective  derivative,  and  by  V  for  the  viscous  derivative,  where: 


v'c-v  =  2(0",  -  o",_1/3) 

The  ADI  sequence  becomes: 

l/+M,ijyA43  =  w?; 

(/  +  -cl)-  hR'-Wtfj] Mtj  =  A$ 


♦n 


«?*  =  fa  +  A,?, 


•4 


(62) 


(63. a) 
(63.6) 
(63. c) 


with  the  time  derivative  taken  at  the  points  placed  Arj/4  from  the  wall  boundary. 
With  the  same  assumption  for  the  r.h.s,  we  get: 


fa  -  -(«&  +  fa  fa  - 


(64) 


6.2.4.  Symmetric  axis  y=0: 

The  partial  derivative  along  the  rj-axis  is  given  by  the  forward  difference  operator  A,  for  the 
convective  derivative: 


*ifaj  =  F",+i  -  fa  (65) 

and  by  A'  for  the  viscous  derivative,  where: 

=  Vffann  -  °U  (66) 

The  vector  solution  qn+i  is  then  given  for  the  symmetric  axis  by  the  equations: 

(Z  +  66.2^1  A?"  =  (67.a) 
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where: 


(68) 


6.2.5.  Inflow  boundary: 

During  all  the  computation,  the  pressure  p ,  the  temperature  T,  and  by  the  state  law  for  a 
perfect  gas,  the  density  p  are  constant  and  equal  to  their  respective  stagnation  values: 

p  =  Po>  T  =  T0,  p  =  p0  (69) 


The  value  of  u  is  updated  at  the  end  of  each  time  step  by  extrapolation  from  the  interior  points 
by  the  mass  conservation  equation  between  the  first  two  sections  of  the  nozzle. 

Then,  the  new  value  of  t;  is  computed  by  the  relation: 

v  =  u  tan  0(y) 

with  the  function  0(y)  given  by  Holder  et  al  [5]. 

Also,  from  the  new  values  of  the  components  u  and  v  of  the  velocity,  the  total  energy  en  is  updated 
by  equation  (4). 


6.3.  Smoothing: 

Finally,  the  numerical  stability  of  this  method  of  simulation  of  a  high  Reynolds  number  flow 
is  here  controlled  by  adding  to  the  r.h.s  of  equation  (49.a)  a  fourth-order  dissipation  term  and  to 
the  l.h.s  of  eqs  (49.a)  and  (49.6)  a  second-order  dissipation  term  [1,8,12].  The  effect  is  for  example 
to  modify  the  ADI  sequence  (55.a,55.6,55.c)  for  the  interior  points  to  the  form: 


[/  +  h6iE?j  -  *,mpJf^V,A,)Jfj]A£;  =  hr?tj  -  J^[(V, ,-A,)2  +  (V,  Ay)*]  ./*«&•  (71.a) 

[/  +  h{8iF?j  -  Re-'S-G^)  -  eimpJr>(VjAj)JiJ}Aq?J  =  A $£.  (71.6) 

C1  =  C;  +  AC,-  (71-c) 

where  the  values  of  the  coefficients  etXp  and  Simp  are  controlled  by  the  linear  stability  condi¬ 
tions.  We  take  here  the  values: 


•'CZp 


—  ®»mp  —  3^* 


(72) 


We  proceed  in  the  same  way  for  the  wall,  the  symmetric  axis  and  the  outflow  boundaries. 
The  operators  A  and  V  are  the  forward  and  backward  difference  operators  defined  before  (see  eqs 


7.  Results: 

The  numerical  method  described  here  has  been  used  in  the  study  of  the  flow  of  a  compressible 
viscous  fluid  in  a  two-dimensional  convergent-divergent  nozzle  (see  figure  1). 

The  numerical  generated  grid  used  for  this  problem  is  shown  in  figure  3.  The  number  of  grid 
points  was  fixed  to  1125  in  order  to  reduce  the  size  of  the  linear  systems  to  be  solved  on  the  parallel 
machine. 

The  mesh  was  of  size  45x25,  where  45  points  were  used  in  the  horizontal  direction  and,  25  points 
in  the  vertical  direction. 

The  combination  of  the  numerical  algorithm,  grid  mapping  and  boundary  conditions  were  here 
as  a  first  step  tested  on  the  VAX  8600,  before  its  application  on  the  iPSC  Intel  Hypercube. 

The  computations  were  stopped  when  the  maximum  relative  change  in  the  Mach  number  in 
the  throat  and  downstream  region  was  sufficiently  small. 

From  the  initial  condition  and  the  inflow  boundary  conditions,  the  Reynolds  number  at  the  be¬ 
ginning  of  this  simulation  was  equal  to  1.5  10®,  being  evaluated  on  the  symmetric  axis,  using  the 
throat  radius  as  the  reference  length. 

The  solution  proposed  was  obtained  for  a  CPU  time  of  about  35  mn.The  main  properties  of 
this  solution  are  here  well  caracterised  by  the  representation  of  the  equi-Mach  number  curbs  along 
the  nozzle  (see  figure  4),  and  by  the  ratio  of  the  pressure  at  the  symmetric  axis  of  each  section  over 
the  stagnation  pressure  (see  figure  5). 

We  get  accurate  results  comparable  to  the  numerical  solution  proposed  by  Cline  in  [3],  In 
addition,  the  capacity  of  adaptability  of  the  numerical  method  used  here  to  a  large  variety  of 
problems  gives  more  power  and  more  interest  to  this  general  approach  [12], 

Also,  we  can  note  that  by  the  definition  of  an  ADI  sequence,  this  numerical  algorithm  has 
intrinsically  a  highly  parallel  character  for  its  application  on  a  distributed  memory  parallel  machine 
[7],  We  propose  now  a  first  reflexion  on  the  parallelization  of  this  problem. 


I 


Figure  5:  Axis  pressure  ratio 


8.  Parallelization: 

Let  us  now  study  the  numerical  simulation  of  this  fluid  flow  on  the  Hypercube,  from  the 
existing  sequential  code. 

In  particular,  we  have  seen  that  the  numerical  algorithm  for  the  solution  of  the  2-D  Navier- 
Stokea  equations  leads  to  an  ADI  sequence,  which  has  already  been  study  both  theorically  and 
pratically  on  the  Hypercube  (see  (7,10|). 

First,  we  recall  the  ADI  sequence  (see  eqs  (71)),  that  we  can  express  here  as  follows: 


AijAW  -  biJ  for  j  —  1, . . 

.,25 

(73.o) 

=  for,' =  2,. 

...45 

(73.6) 

with  the  vector  solution  given  by: 

qntl  =  or*.  +  a<T*. 

(73. c) 

where: 


'  AJj  =  I  +  MtEfj  - 

•  C,"j  =  '  +  WifZ  -  to-'fjCTj)  -  (74) 

for  the  definition  of  EJV,  FA,  GA  of  size  4x4,  and  f(V  see  chapter  VI. 
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The  purpose  of  this  work  is  to  study  the  parallelization  of  the  specific  ADI  sequence  given 
before  (eqs(73))  for  the  computational  domain  of  size  44  x  25  defined  in  chapter  III  (the  first 
column  of  the  initial  domain  of  size  45  x  25,  corresponding  to  the  inflow  boundary  is  computed  by 
extrapolation,  see  chapter  VI). 

More  precisely,  we  solve  the  25  equations  (73.a)  (respectively  44  equations  (73. b))  with  i 
varying  from  2  to  45  (respectively  with  j  varying  from  1  to  25)  which  leads  to  the  computation  of 
25  linear  systems  of  size  176  (44x4)  (respectively  44  linear  systems  of  size  100  (25x4)). 

Pratically,  we  use  a  one  dimensional  domain  decomposition  embedded  on  the  Hypercube  of 
k=2n  processors  where  n  is  the  dimension  of  the  Hypercube. 

The  solution  method  chosen  here  can  be  described  as  follows: 

For  the  first  half  step  of  the  ADI  sequence  (eqs  (73. a)),  the  computational  domain  of  size  44  x  25  is 
decomposed  in  k  horizontal  strips,  where  each  strip  is  assigned  to  one  processor.  Every  processor 
has  then  to  solve  locally  25 /k  linear  systems. 

The  matrix  solution  A gt*”  is  then  transposed  among  the  processors  in  order  to  solve  the  second 
half  step  of  the  ADI  sequence  (eqs  (73. b)),  with  the  maximun  level  of  parallelism,  i.e.  44 /k  systems 
to  solve  per  processor. 

Also,  since  all  linear  systems  are  solved  locally,  and  involve  banded  matrices,  they  can  be 
solved  using  a  band  solver  from  LINPACK. 

The  incremental  solution  Aq?  j  is  moved  back  to  its  original  distribution  for  the  computation 
of  by  equation  (73.c).  Finally,  we  compute  the  updated  values  of  the  matrices  A"^1,  C”tx,  and 
the  vector  .  We  can  then  begin  a  new  time  step  of  the  ADI  sequence  (73.a),  (73.b)  and  (73.c). 

Such  an  approach  has  already  been  made  in  the  study  of  the  Schrodinger  equation  by  an  ADI 
sequence  on  the  iPSC  Intel  Hypercube  (see  [10])  and  will  be  developped  further  for  our  physical 
problem. 

9.  Conclusion: 

An  implicit  finite-difference  serial  computer  program  has  been  developped  to  solve  the  parabolic 
approximation  of  the  2-D  Navier-Stokes  equations  for  a  compressible  fluid. 

Its  application  for  the  simulation  of  a  flow  in  a  convergent-divergent  nozzle  has  provided  an 
accurate  solution.  Further  more,  this  numerical  method,  essentially  based  on  an  ADI  sequence  is 
intrinsically  parallelizable,  which  gives  to  this  general  technique  more  interest  in  the  study  of  fluid 
dynamics  problem,  on  a  parallel  machine  of  the  Hypercube  family. 

The  complete  resolution  of  this  problem  on  the  Hypercube  is  currently  studied  and  will  be 
detailed  in  a  next  paper. 
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